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ABSTRACT 


Kalman filtering techniques are combined with a semiana- 
lytical orbit generator to develop a sequential orbit deter- 
mination algorithm. The algorithm is investigated for com- 
putational efficiency, accuracy, and radius of convergence 
by comparison with truth ephemerides and a Cowell Special 
perturbations filter (GTDS FILTER). Test cases relevant to 
satellite navigation are examined. 

Notation and Symbols 

sub-bar (e.g., x ) = vector 

super-bar (e.g., x) = average or mean; 

also statistical mean 

e (e.g., €?j) = formal indication of the order of the 

quantity 

(e = first, e 2 = second, ...) 

££ = [0 0 0 0 0 1] T 
n = mean motion = yj~r 


Equinoctial Elements 
a = semimajor axis 

h = e sin(o + in) k = e cos(u + IQ) 
p = tan I (l/2) sin Q q = tan 1 (i/2) cos Q 
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X = M + o + If2 = mean longitude 
I = retrograde factor 

super-hat (e.g., x) = predicted estimate 
super-tilde (e.g., x) = updated estimate 
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1. INTRODUCTION 


The current trends in Earth satellite orbit determination 
are toward sequential filtering and onboard computation [ 1 ] . 
The Global Positioning System (GPS) currently employs an 
orbit determination algorithm that updates a batch estimated 
nominal trajectory in real-time with an extended Kalman fil- 
ter [2]. This system is designed to achieve an operational 
accuracy within 1.5 m. Telesat, a satellite communications 
system, has been using a sequential system to support all 
station keeping operations for several years now, with both 
improved accuracy and reduced costs 13]. Many other appli- 
cations exist and will develop for which the timeliness, 
accuracy, and efficiency of a real-time orbit determination 
system are highly desirable. 

Orbit determination processes require two capabilities : 
the ability to accurately propagate an orbit, given initial 
conditions; and some estimation algorithm indicating how 
observations of the satellite should be used in updating the 
ephemeris. Advances in the technology of either capability 
imply corresponding advances in . orbit determination pro- 
cesses. Recently, much work has been done by P. Cefola, et 
al. 14], 15], 16], 17] of CSDL in extending Semianaly tical 
Satellite Theory to allow highly accurate and efficient 
orbit propagation. A. Green 14] developed and used some of 
these results in a batch DC estimation algorithm, finding 
accuracies and convergence properties quite comparable to 
high precision Cowell results. This paper explores the 
implications of these advances in Semianaly tical Satellite 
Theory for sequential orbit determination, considering both 
accuracy and efficiency through comparison with batch and 
sequential filters available from GTDS and Green 14]. 

The organization of the paper is dictated by the struc- 
ture of the orbit determination problem. Summaries of semi- 
analytical satellite theory and sequential filtering are 
presented first. Then their combination into an orbit det- 
ermination algorithm is developed to give the algorithm as 
it was finally implemented. Results are not included here; 
they will be presented at the conference. 


2. SEMI ANALYTICAL SATELLITE THEORY 


The accurate and efficient propagation of an ephemeris 
requires both a precise model of the forces acting on the 
satellite and an accurate and efficient means of integrating 
the equations of motion. The equations of motion are giver, 
by Newton's Second Law as 


d 2 r 

dt 2 


y 

— T 

r 


r = 


^d 


[ 1 ] 
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The terms from left to right are the satellite’s accelera- 
tion, the point-mass gravitational attraction, and all other 
(disturbing) accelerations, due to drag, third bodies, solar 
radiation, etc. The disturbing accelerations are typically 
several orders of magnitude smaller than the point-mass 
force . 

Kou any integrator is most accurate and efficient for 
systems with only small nonlinearities and lou frequencies 
in the force model. Historically, this fact has motivated 
tradeoffs between analytical methods, which use simplified 
force models and analytical approximations to obtain the 
integrated ephemeris efficiently, and numerical methods, 
which retain the full force model and use high precision 
numerical integrators to obtain the integrated ephemeris 
quite accurately. 

To increase the efficiency of an ephemeris generator, it 
is necessary to decrease both the magnitude of the nonli- 
nearities as well as the high frequency content of the force 
model. The magnitude of the nonlinearities can be reduced 
by choice of the orbital elements. For example, Keplerian 
and equinoctial elements incorporate the effects of the 
point mass acceleration, leaving only the disturbing accel- 
eration to be accounted for. The transformation from carte- 
sian position and velocity to such an element set is the 
basis of Gauss' VOP equations. I In the early days of modern 
satellite orbit determination, many element sets incorporat- 
ing different components of the disturbing acceleration were 
experimented with; while they could very efficiently propa- 
gate an ephemeris subject to only their selected perturba- 
tions, to achieve real-world accuracy they had to sacrifice 
all efficiency gains with the inclusion of additional per- 
turbations.] The high frequency content is removed by aver- 
aging these frequencies out; more formally, by transforming 
from the current osculating elements described by the VOP 
equations, to mean elements described by "averaged VOP equa- 
tions." For analytical theories, this whole process was 
done by hand, necessitating simplified force models and 
approximate methods. Semianaly tical satellite theory, 

developed after computers became inexpensive and versatile, 
uses numerical methods to handle those force models that 
cannot be averaged analytically. Since the tradeoff between 
numerical averaging of the force model and the use of a high 
precision integrator on it is in favor of averaging by a 
factor of 10 to 100, semianaly tical satellite theory is much 
more efficient than purely numerical theories. There is one 
problem: the transformation back from the mean elements to 

the osculating elements. The high frequency components or 
short periodics were averaged out and must be recovered 
before the mean elements can be used for anything other than 
long term, approximate prediction. The practical recovery 
of the short periodics constitutes one of the important con- 
tributions of the recent work at CSDL. 
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Semianalytical Satellite Theory at CSDL 


Semianaly tical satellite theory at CSDL is implemented in 
equinoctial elements to avoid singularity problems. The 
basic equations are shown formally in Table I. Key things 
to note are the dependence of the mean element rates on only 
the slowly varying elements a, and the expansion of the 
short periodics ji ( a, T ) as a Fourier series whose coeffi- 
cients similarly depend on only the slowly varying elements 
a. Thus the elements a * and short periodic coefficients 
e Ccr ( a ) and efiaC J[ ) can be and are interpolated, allowing 
efficient evaluation of the osculating elements for many 
output times other than those on the integration grid. This 
is significant since for all averaged theories the computa- 
tional cost is proportional to the number of integration 
steps. Averaging allows large steps, but frequent output 
points could require small steps. 


3. SEBUENTIAL FILTERING THEORY 

The equations of motion for the osculating and mean orbi- 
tal elements are shown in Table I. They are nonlinear, as 
are the equations for range and range rate observations 
given in Table II. The orbit determination problem is to 
estimate the satellite’s orbit given some initial (a priori) 
information and a series of observations over time. It can 
be formulated as an optimal estimation problem: 


estimate x(t) , given 


plant 

observations 


x 


Yv = 


f (x) + w 


X (t ) 

— o 


X 

— o 


h(x(t k ),t k ) + v 


[ 2 ] 


using the y k , such that the variance of the error x - x is 
minimum. x_ Q , w, and v are random and uncorrelated, w and v 
are white noise processes. 

The resulting equations require propagating the probabil- 
ity density function of x Ct) and are very difficult and 
expensive to solve. As a result, most sequential orbit det- 
ermination schemes use some suboptimal filter, usually 
adapted from the Kalman filter, which solves the linear 
optimal estimation problem. The two most common adaptations 
are the Linearised Kalman Filter and the Extended Kalman 
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Table I. The Equations of Semianalytical Satellite Theory* 


Osculating Elements 
Mean elements 

Osculating to mean transformation 
(the near identity transformation) 

Osculating VOP equations 


a* = 

[a,h,k,p,q, X] 

a = 

[a,h,k,p,q] T 

a* = 

[a,h,k,p,q, A] 

a_ 

— «p 

[a,h,k,p,q] 

n 

mi 

a* + er^ (a, A) 

da* 


dt 

n£, + £ F(a, 

— C — — 


Mean VOP equations 


da* _ _ 

-r— = ne_ + e A. (a) 

dt —6 —1 — 


Mean Element Rate 


Short Periodics 


Periodicity of short periodics 


A (a) = f e F (a,A) dX 

-1 - 2t r J q 

e^ilfX) = ^ J [£F (a,X) - e A x (a ) ] dA 


-a 


en u (a,X) ^ dX 


h (a,A + 2lT) = 


r 27T - 

I v* 


A) dX = 0 


Series Expansion of Short 
Periodics 

Assume 


oo 

EF(a,X) = D G 2 0 <a)cos oX 

0=0 _ _ 
+ e z a U)sin aX 


* Extracted from Green [4], which contains a good derivation 
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.27T 


where 


ex a (a) 


= ^ j £F(a,X)dX = EA^a) 


1 r 2ir __ 

£X a (a) = — J £F(£, X) cos OX dX 


zz 


i r 2Tr __ 

Z (a) = - / £F(a r X) sin aX dX 

-a ' tt J Q 

CO 


then 


en (a,X) = Yt ££,(a)sin aX - ee^ ( a) cos oX 
1 a=l 


where 


ec_(a) - 4: z -o ( — } + eD io ( ^ h 

On 2oa 


ESo ( i> - -r e 2 0 ( i> - ec io ( i> % 

On 20a 


Partials 
Solve Vector 


define partials 


State partials equation 


Parameter Partials Equation 


Initial Conditions 


T r“ _T-i 
x = [a* C J 

C = parameter vector in force model 
3 a* (t) 


t(t 'V 




3a* (t ) 
— o 

3a* (t) 


= B„ 


= B„ 


3 c 


_d_ 

dt 


dt 


*(t,t > = [e* — — 1 <Mt, 

° L - * 3a* 3a* J 

'Htrt ) = [e^ — + ^=^1 ¥<t,t 

° L 3l* 3a* J 


t o ) 


3ea, 
+ —1 


$(t,-t ) = i, ¥(t f t ) = o 
o — o — 
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Table II. Range and Range-Rate Satellite Observations 


Orbital elements 

Mean equinoctial elements a* = [a,h,k,p,q,X] T 
Osculating elements a_* = [a,h,k,p,q, X] T 
Cartesian inertial element transformation 


P 


v 


= T (a*) 


Cartesian local tangent element transformation 


= radius to origin of fram on earth's surface 


— LT 


= Dp _r 


= PI + D P 


range 

observation p = /p _ . P TrT , 

—LT — ] LT 


range rate 
observation 



P 

— LT 


v 

-LT 
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Filter. These and other nonlinear filters are discussed in 
Ge lb [ 8 } . 

The Linearized Filter is the most basic adaptation. The 
a priori mean state >T(to) is propagated forward in time to 
generate the nominal trajectory 


V« 




— N o 


= X (t ) 

— o o 


[3] 


The plant and observation equations 12) are then linearized 
about this trajectory to obtain the linear problem 


estimate 

plant 

observation 


Ax(t) , given 

Ax (t) = F (t) Ax (t) + w ; Ax(t Q ) = Ax^ 

Ay k = H(t k ) A * ( V + V 

Ay k « h(x(t k ),t k ) - ^<2^ ( t k )-t k ) 


where 


9f 


F (t ) = 


i» (t) 


[4] 


H< V 


9 h 

W'S: 


The statistics of Ax , u, and v carry over from above. 

—O — 

A Kalman filter can solve the explicit problem l M 1 opti- 
mally, but here the implicit dependence on XN(t) makes the 
solution suboptimal. An Extended Kalman Filter is essen- 
tially a linearized filter that starts over, computing a new 
nominal trajectory, after every observation. Though an 
Extended Filter performs better than a Linearized Filter, 
since the nominal trajectory itself is corrected, the use of 
large step sizes and interpolators for efficiency in the 
semianaly tical ephemeris propagator precludes its use here. 
Rather, a modification of the Linearized Filter will be 
used, as discussed below. The equations for a Linearized 
Kalman Filter are given in Table III. 
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Table III. Linearized Kalman Filter Equations 


Estimation Problem 

x(t) = state to be estimated 
y(t) = scalar observation of x(t) 
w(t) = white state process noise 
v(t) = white observation noise 

plant jc(t) = £(x,t) + w 

observations y(t) = h(x,t) + v at times t 


x(t ) = x 

— o — o 


statistics 


E (w) = 0_ 

E (v) = 0 


E(x ) = x 
— o — o 


E(w(t) W (T) ) = Q6(t “ T) 

E (V (t) V (T) ) = r6 (t - T) 


E (x x ) = P 

— o — o — o 


x , w, v are uncorrelated, 
-o' — ' 


Linearized Kalman Filter Solution 


nominal trajectory 
prediction of estimate and covariance 


x(t) = f(x.t) ; X (t ) = X 

— N -N — N o — o 


transition matrix 

• 

pi 

1 


*<t, Vl ) = 


(x N ,t)| $(t,t._ 1 ) 


•‘VrVi 1 

= I 


state prediction 

A 


A 


Ax(tJ = $(t i ,t i _ 1 ) %(t i _ 1 ) ; 

A 

Ax(t ) = Ax(t. ) 

— 1-1 — 1-1 
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covariance prediction 


P(t.) - ♦(t i ,t i _ 1 > P(t 1 . 1 )* T <t.,t.. 1 ) + A ( t i ft i _ 1 ) 


p(Vi> ■ p, Vi> 

t. 


A(t i ,t i-l ) = / 1 $(t i, T) 2 (T) ^ T(t i' T)dT 


* 1-1 


Update of estimate and covariance 
observation partial 





9 


t.) 

1 


Kalman gain 


observation 


state update 


covariance update 


initialization 


K. = 




H . P (t . ) H . + r 

l li 


Ay (tJ = y(t i ) - h(3^, t ± ) 

A A 

Ax (t . ) = Ax(t.) + K [Ay (t . ) = H Ax (t ) ] 

— X — X — 1 1 11 


P(tj 


(I - K_. H i ) P(tJ 


Ax(t ) = 0 

— o 

'V 

p(t ) = p 

o 
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4 . 


SEMIANALYTICAL KALMAN FILTER DESIGN 


The Kalman Filter equations as given in Table III usually 
allow the means of propagating the nominal trajectory and 
the transition matrices to be arbitrary, since the filter 
only requires the values at observation times. However, 
when optimising the computations for efficiency, the struc- 
tures of the integrator and the filter may become intert- 
wined to produce a more efficient result. This is the case 
for a Semianaly tical Kalman Filter, where the use of inter- 
polators for the state, the transition matrices, and the 
short periodic coefficients has definite implications for 
the overall filter design. 

The Linearised Kalman Filter uses observations over time 
to improve the estimate of a satellite’s orbit. Typically 
the observation times are not known in advance, so the 
underlying ephemeris generator must be able to efficiently 
generate the values of the state and the transition matrices 
at essentially arbitrary times and arbitrarily frequently. 
This requirement does not decrease the efficiency of high 
precision numerical integrators (such as Adams-Cowell , 
etc.), since they are constrained to small step sizes anyway 
and automatically generate the required values at many 
points in time. Analytical and Semianaly tical integrators, 
on the other hand, use very large step sizes, generating the 
required state and transition matrices at only a few points 
in time. Such integrators use interpolators to obtain the 
values at intermediate points in time. The contribution at 
CSDL has been to develop an interpolation method that 
retains the efficiency of analytical integrators and also 
gives values with the accuracies of numerical integrators. 

In the optimization of the Semianaly tical Kalman Filter 
for efficiency, the semianaly tical integrator and the Kalman 
Filter each place requirements on the other. 


The use of interpolators by the integrator over the inte- 
gration grid dictates the use of a Linearized Kalman filter 
inside the integration grid, although the solve vector can 
be updated after processing all the observations in that 
gr id . 

The filter, on the other hand, requires the transition 
matrices $ ( t ^ , t ) , 'Mt£»ti_i ) between adjacent observation 
times t^_^ and t^, . The integrator can most readily supply 
the transition matrices from the beginning of the integra- 
tion grid, <Mti»t 0 )> 'Hti,t 0 ). By using the equations 


♦‘VW ■ ♦‘W ♦'VW 

♦•Wl 1 ■ 

•KVW * ’•'‘W • ♦ ( V t i-i ) ’‘h-i'V 


[ 5 ] 
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we can restate the filter's requirement as supplying 
* (t i»t 0 ). 'Ht i ,^), and $>-l ( ti , t Q ) . While <T 1 (t i ,t 0 ) could 
be calculated directly from <l>(t.i,t 0 ), the expense of comput- 
ing matrix inverses motivates another solution. <f>(t^,t Q ) is 
calculated from a Hermite interpolator using integration 
grid values and rates. Since the rate of 4>“^(t,t 0 ) can be 
calculated as 

$ ^(tjt ) = -4> 1 (t,t ) $(t,t ) $ 1 (t,t ) [6] 

o o o o 


a similar Hermite interpolator can be constructed for 
(t^,t Q ). This interpolator is included in the semiana- 
lytical integrator. 

The last requirement of the filter on the integrator is 
the calculation of A, the contribution of the state process 
noise. Due to the difficulty in defining 2, the process 
noise strength. A, will be calculated as linear in time 


A 


A (t ± - 


t i-l ) 


[7] 


This follows the procedure already incorporated in GTDS [ 9 1 
and appears to work quite well. 

The implementation of the rest of the filter equations is 
straightforward and follows software already in the GTDS 
FILTER subroutines. 

A procedural statement of the final algorithm for imple- 
menting this Semianaly tical Kalman Filter design is given in 
Table IV. 


5. CONCLUSIONS 

An algorithm for implementing a Semianaly tical Kalman 
Filter has been presented. Its implementation is currently 
being completed. Results will be presented at the confer- 
ence . 
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Table IV. The Semi analytical Kalman Filter Algorithm 


Due to use of a Runge kutta integrator , we may consider only one 
integration grid step; all others are processed identically. 


Operations on the Integration Grid 

% * * 

1. t = t update x = x_ + A* 

update P = P 



inititalize Ax = 0 


*(t ,t ) = I 
o o 

V(t f t ) = 0 
o o 

$ _1 (t ft ) = I 
o o 


save in Vs 
save in Os 


2. t = t + At do averaged integration 
o 

A 

obtain x(t), $(t,t ), '{'(tft ) 
— o o 


set up mean interpolators x^, 




3. t = t + At set up interpolators for short periodic 
o 

coefficients EC^ (a; ) , ED^Ca) 
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Operations on the Observation Grid 


1. obtain the new observation, y(t^), at time t = t^. 

A 

2. interpolate for x(t.), $(t.,t ) 'Hit.rt ) 

— 1 1 o l o 

we already have $ ht Q ,t^ ) in 4>s 

3. interpolate for short periodic coefficients 


EC (a (t . ) ) , £D_(a(t.)) 

— o — 1 — O — 1 

4. construct the osculating elements 

_ N __ _ _ _ 

a*(t.) = a*(t.) + Y, sin “ CD (a) cos oA 

1 1 0=1 


5. transform to cartesian elements and construct the nominal observation 

A 

h (x^(t^) ,t^) 
the observation residual 


Ay(tJ = yft.^) - h(x(t i ),t i ) 

and the observation partials 


= 7 <vV -s nb h * B i! b 4 ] 

9x — 


B i = 


3 EH! (a. A) 

3a* 


B 4 = 


3 £]]_ x ( a , A ) 

3 c 


6. Compute the transition matrices 


.-1 


$(t.,t. .) =$(t.,t ) $ (t ,t ) = $(t.,t ) $s 

11-1 lO 1—10 xo 

using $ = (t . n ,t ) , and 4f = ^ (t ,t ) 

S 1-1 O s 1»1 o 
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7. Obtain predicted solve vector and covariance 


Ax(t i ) 




A 


P(t.) 

1 




% 

0 

I 

p( Vi> 


♦'VW 


0 


¥(t i ,t i- 1 ) 


I 


““i'W 


A<t i' t i-1> = A • (t i - w 


8. Complete the update phase of the filter. 
Calculate the gain K. = 


T 

Pit.) H.- 


1 (H. P (t . ) H . + R) 

1 ll 


update the state Ax(t.) = Ax(t.) + K. (Ay(t.) - H. Ax(t.)) 

— 1 — ! 11 1 — 1 

'V * 

update the covariance P(tJ = (I - KH)P(t_^) 


9. Save transition matrices for next observation 

$ = $ 1 (t, f t ) interpolated 

s 1 o 

V = T(t.,t ) interpolated in 2 

s 1 o 


Go to step 1. 
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